#ifndef _CDFLIB_H_
#define _CDFLIB_H_

/*!
 * \file
 * \ingroup pmaths
 * \brief	DCDFLIB Library of Routines for Cumulative Distribution  Functions, Inverses, 
 *			and Other Parameters
 * \details The author accreditation was obtained by searching a UNIX code respositry 
 *			for this source file originally obtained from the PLINK code-base.
 * \authors Barry W. Brown, James Lovato, Kathy Russell 
 * \date	February 1994 
 * \version 1.0
 */

#ifdef __cplusplus
extern "C" {
#endif

/*! \brief Machine precision */
int ipmpar(int *i);

/*! 
* \brief	COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B .GE. 8
* \details	IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY
*			LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X). 
*/
double algdiv(double *a,double *b);

/*!
 * \brief	double precision LN of the GAMMA function
 * \details	<p>If X .le. 6.0, then use recursion to get X below 3
 *			then apply rational approximation number 5236 of
 *			Hart et al, Computer Approximations, John Wiley and
 *			Sons, NY, 1968.</p>
 *			<p>If X .gt. 6.0, then use recursion to get X to at least 12 and
 *		   then use formula 5423 of the same source.</p>
 * \param x value at which scaled log gamma is to be returned. X is DOUBLE PRECISION
 * \returns Returns the natural logarithm of GAMMA(X).
 */
double alngam(double *x);

/*! \brief EVALUATION OF THE FUNCTION LN(1 + A) */
double alnrel(double *a);

/*! 
* \brief	APSER YIELDS THE INCOMPLETE BETA RATIO I(SUB(1-X))(B,A)
* \details	APSER YIELDS THE INCOMPLETE BETA RATIO I(SUB(1-X))(B,A) FOR
*			A .LE. MIN(EPS,EPS*B), B*X .LE. 1, AND X .LE. 0.5. USED WHEN
*			A IS VERY SMALL. USE ONLY IF ABOVE INEQUALITIES ARE SATISFIED.*/
double apser(double *a,double *b,double *x,double *eps);

/*!
* \brief	ASYMPTOTIC EXPANSION FOR IX(A,B) FOR LARGE A AND B.
* \details	LAMBDA = (A + B)*Y - B  AND EPS IS THE TOLERANCE USED.
*			IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT
*			A AND B ARE GREATER THAN OR EQUAL TO 15. */
double basym(double *a,double *b,double *lambda,double *eps);

/*! 
* \brief	EVALUATION OF DEL(A0) + DEL(B0) - DEL(A0 + B0)
* \details	EVALUATION OF  DEL(A0) + DEL(B0) - DEL(A0 + B0)  WHERE
*			LN(GAMMA(A)) = (A - 0.5)*LN(A) - A + 0.5*LN(2*PI) + DEL(A).
*			IT IS ASSUMED THAT A0 .GE. 8 AND B0 .GE. 8. */
double bcorr(double *a0,double *b0);

/*! \brief EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION
 *  \details E = 0.5*LN(2*PI) */
double betaln(double *a0,double *b0);

/*! 
 * \brief	CONTINUED FRACTION EXPANSION FOR IX(A,B) WHEN A,B .GT. 1.
 * \details	IT IS ASSUMED THAT  LAMBDA = (A + B)*Y - B. */
double bfrac(double *a,double *b,double *x,double *y,double *lambda, double *eps);

/*!
\brief ASYMPTOTIC EXPANSION FOR IX(A,B) WHEN A IS LARGER THAN B
\details ASYMPTOTIC EXPANSION FOR IX(A,B) WHEN A IS LARGER THAN B.
     THE RESULT OF THE EXPANSION IS ADDED TO W. IT IS ASSUMED
     THAT A .GE. 15 AND B .LE. 1.  EPS IS THE TOLERANCE USED.
     IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
 */
void bgrat(double *a,double *b,double *x,double *y,double *w, double *eps,int *ierr);

/*! \brief POWER SERIES EXPANSION FOR EVALUATING IX(A,B) 
\details POWER SERIES EXPANSION FOR EVALUATING IX(A,B) WHEN B .LE. 1
     OR B*X .LE. 0.7.  EPS IS THE TOLERANCE USED. */
double bpser(double *a,double *b,double *x,double *eps);

/*! \brief EVALUATION OF THE INCOMPLETE BETA FUNCTION IX(A,B)
\details 
     <p>IT IS ASSUMED THAT A AND B ARE NONNEGATIVE, AND THAT X .LE. 1
     AND Y = 1 - X.  BRATIO ASSIGNS W AND W1 THE VALUES</p>
 
                      <p>W  = IX(A,B)<br>
                      W1 = 1 - IX(A,B)</p>
 
     <p>IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
     IF NO INPUT ERRORS ARE DETECTED THEN IERR IS SET TO 0 AND
     W AND W1 ARE COMPUTED. OTHERWISE, IF AN ERROR IS DETECTED,
     THEN W AND W1 ARE ASSIGNED THE VALUE 0 AND IERR IS SET TO
     ONE OF THE FOLLOWING VALUES ...</p>
 
		<p>
        IERR = 1  IF A OR B IS NEGATIVE<br>
        IERR = 2  IF A = B = 0<br>
        IERR = 3  IF X .LT. 0 OR X .GT. 1<br>
        IERR = 4  IF Y .LT. 0 OR Y .GT. 1<br>
        IERR = 5  IF X + Y .NE. 1<br>
        IERR = 6  IF X = A = 0<br>
        IERR = 7  IF Y = B = 0,</p>
 
\author ALFRED H. MORRIS, JR. (NAVAL SURFACE WARFARE CENTER), DAHLGREN, VIRGINIA
\date  NOV 1991 
*/
void bratio(double *a,double *b,double *x,double *y,double *w, double *w1,int *ierr);

/*! \brief EVALUATION OF  EXP(MU) * (X**A*Y**B/BETA(A,B)) */
double brcmp1(int *mu,double *a,double *b,double *x,double *y);

/*! \brief EVALUATION OF X**A*Y**B/BETA(A,B) */
double brcomp(double *a,double *b,double *x,double *y);

/*! \brief EVALUATION OF IX(A,B) - IX(A+N,B) WHERE N IS A POSITIVE INTEGER. 
\details EPS IS THE TOLERANCE USED. */
double bup(double *a,double *b,double *x,double *y,int *n,double *eps);

/*! \brief Cumulative Distribution Function BETA Distribution
\details <p>Calculates any one parameter of the beta distribution given values for the others.</p>
<p>Cumulative distribution function  (P)  is calculated directly by
     code associated with the following reference.</P>

<p>     DiDinato, A. R. and Morris,  A.   H.  Algorithm 708: Significant
     Digit Computation of the Incomplete  Beta  Function Ratios.  ACM
     Trans. Math.  Softw. 18 (1993), 360-373.</p>

<p>Computation of other parameters involve a seach for a value that
     produces  the desired  value  of P.   The search relies  on  the
     monotinicity of P with the other parameter.</p>
<p> The beta density is proportional to t^(A-1) * (1-t)^(B-1)>/p>
\param which Integer indicating which of the next four argument<br>
               values is to be calculated from the others.<br>
               Legal range: 1..4<br>
               iwhich = 1 : Calculate P and Q from X,Y,A and B<br>
               iwhich = 2 : Calculate X and Y from P,Q,A and B<br>
               iwhich = 3 : Calculate A from P,Q,X,Y and B<br>
               iwhich = 4 : Calculate B from P,Q,X,Y and A<br>
\param p The integral from 0 to X of the chi-square distribution. Input range: [0, 1].
\param q 1-P. Input range: [0, 1]. P + Q = 1.0.
\param x Upper limit of integration of beta density. Input range: [0,1]. Search range: [0,1]
\param y 1-X. Input range: [0,1]. Search range: [0,1], X + Y = 1.0.
\param a The first parameter of the beta density.  Input range: (0, +infinity). Search range: [1D-300,1D300]
\param b The second parameter of the beta density.
            Input range: (0, +infinity).
            Search range: [1D-300,1D300]
\param status 0 if calculation completed correctly<br>
               -I if input parameter number I is out of range<br>
                1 if answer appears to be lower than lowest search bound<br>
                2 if answer appears to be higher than greatest
                  search bound<br>
                3 if P + Q .ne. 1<br>
                4 if X + Y .ne. 1<br>
\param bound Undefined if STATUS is 0.
Bound exceeded by parameter number I if STATUS is negative.

               Lower search bound if STATUS is 1.
               Upper search bound if STATUS is 2.
*/
void cdfbet(int *which,double *p,double *q,double *x,double *y, double *a,double *b,int *status,double *bound);

/*!
\brief Cumulative Distribution Function BINomial distribution
\details
     <p>Calculates any one parameter of the binomial
     distribution given values for the others.</p>

	 <p> Formula  26.5.24    of   Abramowitz  and    Stegun,  Handbook   of
     Mathematical   Functions (1966) is   used  to reduce the  binomial
     distribution  to  the  cumulative incomplete    beta distribution.</p>

     <P>Computation of other parameters involve a seach for a value that
     produces  the desired  value  of P.   The search relies  on  the
     monotinicity of P with the other parameter.</P>

\param which                              
			   Integer indicating which of the next four argument
               values is to be calculated from the others.<br>
               Legal range: 1..4<br>
               iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR<br>
               iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR<br>
               iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR<br>
               iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN<br>
\param p
		    The cumulation from 0 to S of the binomial distribution.
            (Probablility of S or fewer successes in XN trials each
            with probability of success PR.)<br>
            Input range: [0,1].
\param q
            1-P.<br>
            Input range: [0, 1].<br>
            P + Q = 1.0.<br>
\param s
			The number of successes observed.<br>
            Input range: [0, XN]<br>
            Search range: [0, XN]<br>
\param xn 
			  The number of binomial trials.<br>
              Input range: (0, +infinity).<br>
              Search range: [1E-300, 1E300]<br>
\param pr
              The probability of success in each binomial trial.<br>
              Input range: [0,1].<br>
              Search range: [0,1]<br>
\param ompr 1-PR
              Input range: [0,1].<br>
              Search range: [0,1]<br>
              PR + OMPR = 1.0<br>
\param status Status Flag<br>
                0 if calculation completed correctly<br>
               -I if input parameter number I is out of range<br>
                1 if answer appears to be lower than lowest
                  search bound<br>
                2 if answer appears to be higher than greatest
                  search bound<br>
                3 if P + Q .ne. 1<br>
                4 if PR + OMPR .ne. 1<br>
\param bound Undefined if STATUS is 0<br>

               Bound exceeded by parameter number I if STATUS
               is negative.<br>
               Lower search bound if STATUS is 1.<br>
               Upper search bound if STATUS is 2.<br> */
void cdfbin(int *which,double *p,double *q,double *s,double *xn, double *pr,double *ompr,int *status,double *bound);

/*!
\brief Cumulative Distribution Function CHI-Square distribution
\details 
     <p>Calculates any one parameter of the chi-square
     distribution given values for the others.</p>

	 <p>Formula    26.4.19   of Abramowitz  and     Stegun, Handbook  of
     Mathematical Functions   (1966) is used   to reduce the chisqure
     distribution to the incomplete distribution.</p>

     <p>Computation of other parameters involve a seach for a value that
     produces  the desired  value  of P.   The search relies  on  the
     monotinicity of P with the other parameter.</p>

\param which
	           Integer indicating which of the next three argument
               values is to be calculated from the others.<br>
               Legal range: 1..3<br>
               iwhich = 1 : Calculate P and Q from X and DF<br>
               iwhich = 2 : Calculate X from P,Q and DF<br>
               iwhich = 3 : Calculate DF from P,Q and X<br>
\param p
            The integral from 0 to X of the chi-square
            distribution.<br>
            Input range: [0, 1].
\param q
            1-P.<br>
            Input range: (0, 1].<br>
            P + Q = 1.0.<br>
\param x
            Upper limit of integration of the non-central
            chi-square distribution.<br>
            Input range: [0, +infinity).<br>
            Search range: [0,1E300]
\param df 
             Degrees of freedom of the
             chi-square distribution.<br>
             Input range: (0, +infinity).<br>
             Search range: [ 1E-300, 1E300]

\param status 
				Status Flag<br>
                0 if calculation completed correctly<br>
               -I if input parameter number I is out of range<br>
                1 if answer appears to be lower than lowest
                  search bound<br>
                2 if answer appears to be higher than greatest
                  search bound<br>
                3 if P + Q .ne. 1<br>
               10 indicates error returned from cumgam.  See
                  references in cdfgam<br>
\param bound
			   Undefined if STATUS is 0
               Bound exceeded by parameter number I if STATUS
               is negative.
               Lower search bound if STATUS is 1.
               Upper search bound if STATUS is 2. */
void cdfchi(int *which,double *p,double *q,double *x,double *df, int *status,double *bound);

/*!
\brief Cumulative Distribution Function Non-central Chi-Square
\details     
     <p>Calculates any one parameter of the non-central chi-square
     distribution given values for the others.</p>

	 <p>Formula  26.4.25   of   Abramowitz   and   Stegun,  Handbook  of
     Mathematical  Functions (1966) is used to compute the cumulative
     distribution function.</p>

     <p>Computation of other parameters involve a seach for a value that
     produces  the desired  value  of P.   The search relies  on  the
     monotinicity of P with the other parameter.</p>

	 <p><STRONG>Warning:-</strong><br>
	 The computation time  required for this  routine is proportional
     to the noncentrality  parameter  (PNONC).  Very large  values of
     this parameter can consume immense  computer resources.  This is
     why the search range is bounded by 10,000.</p>

\param which
		Integer indicating which of the next three argument
               values is to be calculated from the others.<br>
               Input range: 1..4<br>
               iwhich = 1 : Calculate P and Q from X and DF<br>
               iwhich = 2 : Calculate X from P,DF and PNONC<br>
               iwhich = 3 : Calculate DF from P,X and PNONC<br>
               iwhich = 3 : Calculate PNONC from P,X and DF<br>
\param p
			The integral from 0 to X of the non-central chi-square
            distribution.<br>
            Input range: [0, 1-1E-16).<br>
\param q
			1-P.
            Q is not used by this subroutine and is only included
            for similarity with other cdf* routines.
\param x
		    Upper limit of integration of the non-central
            chi-square distribution.<br>
            Input range: [0, +infinity).<br>
            Search range: [0,1E300]
\param df
             Degrees of freedom of the non-central
             chi-square distribution.<br>
             Input range: (0, +infinity).<br>
             Search range: [ 1E-300, 1E300]
\param pnonc
                Non-centrality parameter of the non-central
                chi-square distribution.<br>
                Input range: [0, +infinity).<br>
                Search range: [0,1E4]
\param status Status Flag
				0 if calculation completed correctly<br>
               -I if input parameter number I is out of range<br>
                1 if answer appears to be lower than lowest
                  search bound<br>
                2 if answer appears to be higher than greatest
                  search bound<br> 
\param bound   Undefined if STATUS is 0<br>
               Bound exceeded by parameter number I if STATUS
               is negative.<br>
               Lower search bound if STATUS is 1.<br>
               Upper search bound if STATUS is 2.<br>
*/
void cdfchn(int *which,double *p,double *q,double *x,double *df, double *pnonc,int *status,double *bound);

/*!
\brief Cumulative Distribution Function F distribution
\details
     <p>Calculates any one parameter of the F distribution
     given values for the others.</p>

	 <p>Formula   26.6.2   of   Abramowitz   and   Stegun,  Handbook  of
     Mathematical  Functions (1966) is used to reduce the computation
     of the  cumulative  distribution function for the  F  variate to
     that of an incomplete beta.</p>

     <p>Computation of other parameters involve a seach for a value that
     produces  the desired  value  of P.   The search relies  on  the
     monotinicity of P with the other parameter.</p>
        
	 <p><strong>WARNING:-</strong><br>
	 The value of the  cumulative  F distribution is  not necessarily
     monotone in  either degrees of freedom.  There  thus may  be two
     values  that  provide a given CDF  value.   This routine assumes
     monotonicity and will find an arbitrary one of the two values.</p>

\param which
			   Integer indicating which of the next four argument
               values is to be calculated from the others.
               Legal range: 1..4
               iwhich = 1 : Calculate P and Q from F,DFN and DFD
               iwhich = 2 : Calculate F from P,Q,DFN and DFD
               iwhich = 3 : Calculate DFN from P,Q,F and DFD
               iwhich = 4 : Calculate DFD from P,Q,F and DFN
\param p
        The integral from 0 to F of the f-density.<br>
              Input range: [0,1].
\param  q
        1-P.<br>
              Input range: (0, 1].<br>
              P + Q = 1.0.
\param  f
        Upper limit of integration of the f-density.<br>
              Input range: [0, +infinity).<br>
              Search range: [0,1E300]
\param dfn
      Degrees of freedom of the numerator sum of squares.<br>
               Input range: (0, +infinity).<br>
               Search range: [ 1E-300, 1E300]
\param dfd
     Degrees of freedom of the denominator sum of squares.<br>
               Input range: (0, +infinity).<br>
               Search range: [ 1E-300, 1E300]
\param status Status Flag<br>
				0 if calculation completed correctly<br>
               -I if input parameter number I is out of range<br>
                1 if answer appears to be lower than lowest
                  search bound<br>
                2 if answer appears to be higher than greatest
                  search bound<br>
                3 if P + Q .ne. 1<br>
\param bound
			   Undefined if STATUS is 0
               Bound exceeded by parameter number I if STATUS
               is negative.
               Lower search bound if STATUS is 1.
               Upper search bound if STATUS is 2. */
void cdff(int *which,double *p,double *q,double *f,double *dfn, double *dfd,int *status,double *bound);

/*! 
\brief Cumulative Distribution Function Non-central F distribution
\details
	<p>Calculates any one parameter of the Non-central F
    distribution given values for the others.</p>

	<p>Formula  26.6.20   of   Abramowitz   and   Stegun,  Handbook  of
     Mathematical  Functions (1966) is used to compute the cumulative
     distribution function.</p>

     <p>Computation of other parameters involve a seach for a value that
     produces  the desired  value  of P.   The search relies  on  the
     monotinicity of P with the other parameter.</p>

     <p>
	 <strong>WARNING:-</strong><br>
	 The computation time  required for this  routine is proportional
     to the noncentrality  parameter  (PNONC).  Very large  values of
     this parameter can consume immense  computer resources.  This is
     why the search range is bounded by 10,000.</p>

     <p>
	 <strong>WARNING:-</strong><br>
	 The  value  of the  cumulative  noncentral F distribution is not
     necessarily monotone in either degrees  of freedom.  There  thus
     may be two values that provide a given  CDF value.  This routine
     assumes monotonicity  and will find  an arbitrary one of the two</p>
   
\param
     which --> Integer indicating which of the next five argument
               values is to be calculated from the others.<br>
               Legal range: 1..5<br>
               iwhich = 1 : Calculate P and Q from F,DFN,DFD and PNONC<br>
               iwhich = 2 : Calculate F from P,Q,DFN,DFD and PNONC<br>
               iwhich = 3 : Calculate DFN from P,Q,F,DFD and PNONC<br>
               iwhich = 4 : Calculate DFD from P,Q,F,DFN and PNONC<br>
               iwhich = 5 : Calculate PNONC from P,Q,F,DFN and DFD<br>

\param p
       The integral from 0 to F of the non-central f-density.<br>
              Input range: [0,1-1E-16).
\param q
       1-P.<br>
              Q is not used by this subroutine and is only included
              for similarity with other cdf* routines.
\param f
       Upper limit of integration of the non-central f-density.<br>
              Input range: [0, +infinity).<br>
              Search range: [0,1E300]
\param dfn
     Degrees of freedom of the numerator sum of squares.<br>
               Input range: (0, +infinity).<br>
               Search range: [ 1E-300, 1E300]
\param dfd
     Degrees of freedom of the denominator sum of squares.<br>
               Must be in range: (0, +infinity).<br>
               Input range: (0, +infinity).<br>
               Search range: [ 1E-300, 1E300]
\param phonc
     The non-centrality parameter<br>
               Input range: [0,infinity)<br>
               Search range: [0,1E4]
\param status Status flag<br>
				0 if calculation completed correctly<br>
               -I if input parameter number I is out of range<br>
                1 if answer appears to be lower than lowest
                  search bound<br>
                2 if answer appears to be higher than greatest
                  search bound<br>
                3 if P + Q .ne. 1<br>
\param bound
			   Undefined if STATUS is 0
               Bound exceeded by parameter number I if STATUS
               is negative.
               Lower search bound if STATUS is 1.
               Upper search bound if STATUS is 2. */
void cdffnc(int *which,double *p,double *q,double *f,double *dfn, double *dfd,double *phonc,int *status,double *bound);

/*!
\brief Cumulative Distribution Function GAMma Distribution
\details
     <p>Calculates any one parameter of the gamma
     distribution given values for the others.</p>

	 <p>Cumulative distribution function (P) is calculated directly by
     the code associated with:</p>

     <p>DiDinato, A. R. and Morris, A. H. Computation of the  incomplete
     gamma function  ratios  and their  inverse.   ACM  Trans.  Math.
     Softw. 12 (1986), 377-393.</p>

     <p>Computation of other parameters involve a seach for a value that
     produces  the desired  value  of P.   The search relies  on  the
     monotinicity of P with the other parameter.</p>

	 <p>The gamma density is proportional to
       T**(SHAPE - 1) * EXP(- SCALE * T)</p>

\param which --> Integer indicating which of the next four argument
               values is to be calculated from the others.<br>
               Legal range: 1..4<br>
               iwhich = 1 : Calculate P and Q from X,SHAPE and SCALE<br>
               iwhich = 2 : Calculate X from P,Q,SHAPE and SCALE<br>
               iwhich = 3 : Calculate SHAPE from P,Q,X and SCALE<br>
               iwhich = 4 : Calculate SCALE from P,Q,X and SHAPE<br>

\param p
            The integral from 0 to X of the gamma density.<br>
            Input range: [0,1].
\param q
            Input range: (0, 1].<br>
            P + Q = 1.0.
\param x
     The upper limit of integration of the gamma density.<br>
            Input range: [0, +infinity).<br>
            Search range: [0,1E300]
\param shape
     The shape parameter of the gamma density.<br>
                Input range: (0, +infinity).<br>
                Search range: [1E-300,1E300]
\param scale
     The scale parameter of the gamma density.<br>
                Input range: (0, +infinity).<br>
                Search range: (1E-300,1E300]
\param status Status Flag<br>
				0 if calculation completed correctly<br>
               -I if input parameter number I is out of range<br>
                1 if answer appears to be lower than lowest
                  search bound<br>
                2 if answer appears to be higher than greatest
                  search bound<br>
                3 if P + Q .ne. 1<br>
                10 if the gamma or inverse gamma routine cannot
                   compute the answer.  Usually happens only for
                   X and SHAPE very large (gt 1E10 or more)<br>
\param bound
			   Undefined if STATUS is 0
               Bound exceeded by parameter number I if STATUS
               is negative.
               Lower search bound if STATUS is 1.
               Upper search bound if STATUS is 2.
*/
void cdfgam(int *which,double *p,double *q,double *x,double *shape, double *scale,int *status,double *bound);

/*!
\brief Cumulative Distribution Function Negative BiNomial distribution
\details
     <p>Calculates any one parameter of the negative binomial
     distribution given values for the others.</p>

     <p>The  cumulative  negative   binomial  distribution  returns  the
     probability that there  will be  F or fewer failures before  the
     XNth success in binomial trials each of which has probability of
     success PR.<p>

     <p>The individual term of the negative binomial is the probability of
     S failures before XN successes and is
          Choose( S, XN+S-1 ) * PR^(XN) * (1-PR)^S</p>

     <p>Formula   26.5.26   of   Abramowitz  and  Stegun,  Handbook   of
     Mathematical Functions (1966) is used  to  reduce calculation of
     the cumulative distribution  function to that of  an  incomplete
     beta.</p>

     <p>Computation of other parameters involve a seach for a value that
     produces  the desired  value  of P.   The search relies  on  the
     monotinicity of P with the other parameter.</p>

\param which
     Integer indicating which of the next four argument
               values is to be calculated from the others.<br>
               Legal range: 1..4<br>
               iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR<br>
               iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR<br>
               iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR<br>
               iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN<br>
\param p
		    The cumulation from 0 to S of the  negative
            binomial distribution.<br>
            Input range: [0,1].
\param q
			1-P.<br>
            Input range: (0, 1].<br>
            P + Q = 1.0.
\param s
			The upper limit of cumulation of the binomial distribution.<br>
            There are F or fewer failures before the XNth success.<br>
            Input range: [0, +infinity).<br>
            Search range: [0, 1E300]
\param xn
			  The number of successes.<br>
              Input range: [0, +infinity).<br>
              Search range: [0, 1E300]
\param pr
			  The probability of success in each binomial trial.<br>
              Input range: [0,1].<br>
              Search range: [0,1].
\param ompr
			  1-PR<br>
              Input range: [0,1].<br>
              Search range: [0,1]<br>
              PR + OMPR = 1.0
\param status Status Flag
				0 if calculation completed correctly<br>
               -I if input parameter number I is out of range<br>
                1 if answer appears to be lower than lowest
                  search bound<br>
                2 if answer appears to be higher than greatest
                  search bound<br>
                3 if P + Q .ne. 1<br>
                4 if PR + OMPR .ne. 1<br>
\param bound
               Undefined if STATUS is 0
               Bound exceeded by parameter number I if STATUS
               is negative.
               Lower search bound if STATUS is 1.
               Upper search bound if STATUS is 2. */
void cdfnbn(int *which,double *p,double *q,double *s,double *xn, double *pr,double *ompr,int *status,double *bound);

/*!
\brief Cumulative Distribution Function NORmal distribution
     <p>Calculates any one parameter of the normal
     distribution given values for the others.</p>

	 <p>A slightly modified version of ANORM from:-</p>

     <p>Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
     Package of Special Function Routines and Test Drivers"
     acm Transactions on Mathematical Software. 19, 22-32.</p>

     <p>is used to calulate the  cumulative standard normal distribution.</p>

     <p>The rational functions from pages  90-95  of Kennedy and Gentle,
     Statistical  Computing,  Marcel  Dekker, NY,  1980 are  used  as
     starting values to Newton's Iterations which compute the inverse
     standard normal.  Therefore no  searches  are necessary for  any
     parameter.</P>

     <p>For X < -15, the asymptotic expansion for the normal is used  as
     the starting value in finding the inverse standard normal.
     This is formula 26.2.12 of Abramowitz and Stegun.</p>

	 <p>N/B The normal density is proportional to
      exp( - 0.5 * (( X - MEAN)/SD)**2)</p>

\param which  Integer indicating  which of the  next  parameter
     values is to be calculated using values  of the others.
	           <ul>
			   <li>Legal range: 1..4</li>
               <li>which = 1 : Calculate P and Q from X,MEAN and SD</li>
               <li>which = 2 : Calculate X from P,Q,MEAN and SD</li>
               <li>which = 3 : Calculate MEAN from P,Q,X and SD</li>
               <li>which = 4 : Calculate SD from P,Q,X and MEAN</li>
			   </ul>
\param p 
            The integral from -infinity to X of the normal density.
            Input range: (0,1].
\param q
            1-P.
            Input range: (0, 1].
            P + Q = 1.0.
\param x
             Upper limit of integration of the normal-density.
             Input range: ( -infinity, +infinity)
\param mean
               The mean of the normal density.
               Input range: (-infinity, +infinity)
\param sd
             Standard Deviation of the normal density.
             Input range: (0, +infinity).
\param status Status Flag
				<ul>
				<li>0 if calculation completed correctly</li>
               <li>-I if input parameter number I is out of range</li>
                <li>1 if answer appears to be lower than lowest
                  search bound</li>
                <li>2 if answer appears to be higher than greatest
                  search bound</li>
                <li>3 if P + Q .ne. 1</li>
				</ul>
\param bound
               Undefined if STATUS is 0
               Bound exceeded by parameter number I if STATUS
               is negative.
               Lower search bound if STATUS is 1.
               Upper search bound if STATUS is 2.
*/
void cdfnor(int *which,double *p,double *q,double *x,double *mean, double *sd,int *status,double *bound);

/*!
\brief cumulative Distribution Function POIsson distribution
\details     
     <P>Calculates any one parameter of the Poisson
     distribution given values for the others.</P>
     <p>Formula   26.4.21  of   Abramowitz  and   Stegun,   Handbook  of
     Mathematical Functions (1966) is used  to reduce the computation
     of  the cumulative distribution function to that  of computing a
     chi-square, hence an incomplete gamma function.</p>

     <p>Cumulative  distribution function  (P) is  calculated  directly.
     Computation of other parameters involve a seach for a value that
     produces  the desired value of  P.   The  search relies  on  the
     monotinicity of P with the other parameter.</p>

\param which Integer indicating which  argument
               value is to be calculated from the others.
			   <ul>
               <li>Legal range: 1..3</li>
               <li>iwhich = 1 : Calculate P and Q from S and XLAM</li>
               <li>iwhich = 2 : Calculate A from P,Q and XLAM</li>
               <li>iwhich = 3 : Calculate XLAM from P,Q and S</li>
			   </ul>
\param p The cumulation from 0 to S of the poisson density.
               Input range: [0,1].
\param q 1-P.
               Input range: (0, 1].
               P + Q = 1.0.

\param s       Upper limit of cumulation of the Poisson.
               Input range: [0, +infinity).
               Search range: [0,1E300]

\param xlam    Mean of the Poisson distribution.
               Input range: [0, +infinity).
               Search range: [0,1E300]

\param status Status Flags 
				<ul>
				<li>0 if calculation completed correctly</li>
               <li>-I if input parameter number I is out of range</li>
                <li>1 if answer appears to be lower than lowest
                  search bound</li>
                <li>2 if answer appears to be higher than greatest
                  search bound</li>
                <li>3 if P + Q .ne. 1</li>
				</ul>

\param bound   Undefined if STATUS is 0.
               Bound exceeded by parameter number I if STATUS
               is negative.
               Lower search bound if STATUS is 1.
               Upper search bound if STATUS is 2. */
void cdfpoi(int *which,double *p,double *q,double *s,double *xlam, int *status,double *bound);


/*! 
\brief Cumulative Distribution Function T distribution
\details     
     <p>Calculates any one parameter of the t distribution given
     values for the others.</p>
     
	 <p>Formula  26.5.27  of   Abramowitz   and  Stegun,   Handbook   of
     Mathematical Functions  (1966) is used to reduce the computation
     of the cumulative distribution function to that of an incomplete
     beta.</p>

     <p>Computation of other parameters involve a seach for a value that
     produces  the desired  value  of P.   The search relies  on  the
     monotinicity of P with the other parameter.</p>

\param which Integer indicating which  argument
               values is to be calculated from the others.<br>
               Legal range: 1..3<br>
               iwhich = 1 : Calculate P and Q from T and DF<br>
               iwhich = 2 : Calculate T from P,Q and DF<br>
               iwhich = 3 : Calculate DF from P,Q and T
\param p The integral from -infinity to t of the t-density.
               Input range: (0,1].
\param q 1-P. <br> Input range: (0, 1].
               P + Q = 1.0.
\param t Upper limit of integration of the t-density.<br>
         Input range: ( -infinity, +infinity).<br>
         Search range: [ -1E300, 1E300 ]
\param df Degrees of freedom of the t-distribution.<br>
                Input range: (0 , +infinity).<br>
                Search range: [1e-300, 1E10]
\param status Status flag<br>
                0 if calculation completed correctly<br>
               -I if input parameter number I is out of range<br>
                1 if answer appears to be lower than lowest
                  search bound<br>
                2 if answer appears to be higher than greatest
                  search bound<br>
                3 if P + Q .ne. 1<br>
\param bound Undefined if STATUS is 0<br>

               Bound exceeded by parameter number I if STATUS
               is negative.<br>

               Lower search bound if STATUS is 1.<br>

               Upper search bound if STATUS is 2.<br> */
void cdft(int *which,double *p,double *q,double *t,double *df, int *status,double *bound);

/*! \brief Double precision cUMulative incomplete BETa distribution
 \details
     <p>Calculates the cdf to X of the incomplete beta distribution
     with parameters a and b.  This is the integral from 0 to x
     of (1/B(a,b))*f(t)) where f(t) = t**(a-1) * (1-t)**(b-1)</p>
	 <p>Calls the routine bratio().</p>
 
\param x Upper limit of integration.
\param y 1 - x                             
\param a First parameter of the beta distribution.
\param b Second parameter of the beta distribution.
\param cum Cumulative incomplete beta distribution.
\param ccum Compliment of Cumulative incomplete beta distribution.                                      
\author 
     Didonato, Armido R. and Morris, Alfred H. Jr. (1992) Algorithim
     708 Significant Digit Computation of the Incomplete Beta Function
     Ratios. ACM ToMS, Vol.18, No. 3, Sept. 1992, 360-373. */
void cumbet(double *x,double *y,double *a,double *b,double *cum,double *ccum);

/*! \brief CUmulative BINomial distribution
\details
     Returns the probability   of 0  to  S  successes in  XN   binomial
     trials, each of which has a probability of success, PBIN.
\param s The upper limit of cumulation of the binomial distribution.
\param xn The number of binomial trials.
\param pr The probability of success in each binomial trial.
\param ompr 1 - PBIN
\param cum Cumulative binomial distribution.
\param ccum Compliment of Cumulative binomial distribution. 
\author
	Formula  26.5.24    of   Abramowitz  and    Stegun,  Handbook   of
     Mathematical   Functions (1966) is   used  to reduce the  binomial
     distribution  to  the  cumulative    beta distribution. */
void cumbin(double *s,double *xn,double *pr,double *ompr, double *cum,double *ccum);

/*! \brief CUMulative of the CHi-square distribution
\details 
     Calculates the cumulative chi-square distribution.
	 Calls incomplete gamma function cumgam().
\param x Upper limit of integration of the
                 chi-square distribution.
\param df Degrees of freedom of the
                 chi-square distribution.
\param cum Cumulative chi-square distribution.
\param ccum Compliment of Cumulative chi-square distribution. */
void cumchi(double *x,double *df,double *cum,double *ccum);

/*!
\brief CUMulative of the Non-central CHi-square distribution
 \details
     <p>Calculates     the       cumulative      non-central    chi-square
     distribution, i.e.,  the probability   that  a   random   variable
     which    follows  the  non-central chi-square  distribution,  with
     non-centrality  parameter    PNONC  and   continuous  degrees   of
     freedom DF, is less than or equal to X.</p>
	 <p>Local variables:-<br>
	 <pre>
	 EPS     --- Convergence criterion.  The sum stops when a
                 term is less than EPS*SUM.
                                                 EPS is DOUBLE PRECISIO
 
     NTIRED  --- Maximum number of terms to be evaluated
                 in each sum.
                                                 NTIRED is INTEGER
 
     QCONV   --- .TRUE. if convergence achieved -
                 i.e., program did not stop on NTIRED criterion.
                                                 QCONV is LOGICAL
 
     CCUM <-- Compliment of Cumulative non-central
              chi-square distribution.
	 </pre>
	 </p>
\param x Upper limit of integration of the non-central
                 chi-square distribution.
\param df Degrees of freedom of the non-central
                 chi-square distribution.
\param pnonc Non-centrality parameter of the non-central
                 chi-square distribution.
\param cum Cumulative non-central chi-square distribution.
\param ccum Compliment of Cumulative non-central chi-square distribution

\author
	 Uses  formula  26.4.25   of  Abramowitz  and  Stegun, Handbook  of
     Mathematical    Functions,  US   NBS   (1966)    to calculate  the
     non-central chi-square. */
void cumchn(double *x,double *df,double *pnonc,double *cum, double *ccum);

/*! \brief CUMulative F distribution
 \details 
     <p>Computes  the  integral from  0  to  F of  the f-density  with DFN
     and DFD degrees of freedom.</p>
	<p>N/B If F is less than or equal to 0, 0 is returned.</p>
 
\param f Upper limit of integration of the f-density.
\param dfn Degrees of freedom of the numerator sum of squares.
\param dfd Degrees of freedom of the denominator sum of squares.
\param cum Cumulative f distribution.
\param ccum Compliment of Cumulative f distribution.
\author 
     Formula  26.5.28 of  Abramowitz and   Stegun   is  used to  reduce
     the cumulative F to a cumulative beta distribution. */
void cumf(double *f,double *dfn,double *dfd,double *cum,double *ccum);

/*! \brief F -NON- -C-ENTRAL F DISTRIBUTION
\details
     <p>COMPUTES NONCENTRAL F DISTRIBUTION WITH DFN AND DFD
     DEGREES OF FREEDOM AND NONCENTRALITY PARAMETER PNONC</p>

	 <p>USES FORMULA 26.6.20 OF REFERENCE FOR INFINITE SERIES.
     SERIES IS CALCULATED BACKWARD AND FORWARD FROM J = LAMBDA/2
     (THIS IS THE TERM WITH THE LARGEST POISSON WEIGHT) UNTIL
     THE CONVERGENCE CRITERION IS MET.</p>
 
     <p>FOR SPEED, THE INCOMPLETE BETA FUNCTIONS ARE EVALUATED
     BY FORMULA 26.5.16.</p>

	 <p>N/B THE SUM CONTINUES UNTIL A SUCCEEDING TERM IS LESS THAN EPS
     TIMES THE SUM (OR THE SUM IS LESS THAN 1.0E-20).  EPS IS
     SET TO 1.0E-4 IN A DATA STATEMENT WHICH CAN BE CHANGED.</p>
 
\param f UPPER LIMIT OF INTEGRATION OF NONCENTRAL F IN EQUATION
\param dfn DEGREES OF FREEDOM OF NUMERATOR
\param dfd DEGREES OF FREEDOM OF DENOMINATOR
\param pnonc NONCENTRALITY PARAMETER.
\param cum CUMULATIVE NONCENTRAL F DISTRIBUTION
\param ccum COMPLIMENT OF CUMMULATIVE
\author 
 <pre>
     HANDBOOD OF MATHEMATICAL FUNCTIONS
     EDITED BY MILTON ABRAMOWITZ AND IRENE A. STEGUN
     NATIONAL BUREAU OF STANDARDS APPLIED MATEMATICS SERIES - 55
     MARCH 1965
     P 947, EQUATIONS 26.6.17, 26.6.18</pre> */
void cumfnc(double *f,double *dfn,double *dfd,double *pnonc, double *cum,double *ccum);

 /* !\brief Double precision cUMulative incomplete GAMma distribution
 \details    
     <p>Computes   the  cumulative        of    the     incomplete   gamma
     distribution, i.e., the integral from 0 to X of
          (1/GAM(A))*EXP(-T)*T**(A-1) DT
     where GAM(A) is the complete gamma function of A, i.e.,
          GAM(A) = integral from 0 to infinity of
                    EXP(-T)*T**(A-1) DT</p>
<p>Calls the routine gratio()</p>
 \param x The upper limit of integration of the incomplete gamma.
 \param a The shape parameter of the incomplete gamma.
 \param cum Cumulative incomplete gamma distribution.
 \param ccum Compliment of Cumulative incomplete gamma distribution.*/
 void cumgam(double *x,double *a,double *cum,double *ccum);

/*! \brief CUmulative Negative BINomial distribution
\details 
     <p>Returns the probability that it there will be S or fewer failures
     before there are XN successes, with each binomial trial having
     a probability of success PR.<p>
 
 <p><pre>
     Prob(# failures = S | XN successes, PR)  =
                        ( XN + S - 1 )
                        (            ) * PR^XN * (1-PR)^S
                        (      S     )
 </pre></p>
 \param s The number of failures
 \param xn The number of successes
 \param pr The probability of success in each binomial trial.
 \param ompr 1 - PR
 \param cum Cumulative negative binomial distribution.
 \param ccum Compliment of Cumulative negative binomial distribution.
 \author                                                  
     Formula  26.5.26    of   Abramowitz  and    Stegun,  Handbook   of
     Mathematical   Functions (1966) is   used  to reduce the  negative
     binomial distribution to the cumulative beta distribution. */
void cumnbn(double *s,double *xn,double *pr,double *ompr, double *cum,double *ccum);

/*! \brief Computes the cumulative  of    the  normal   distribution
\details
<p>Computes the cumulative  of    the  normal   distribution,   i.e.,
     the integral from -infinity to x of
          (1/sqrt(2*pi)) exp(-u*u/2) du</p>
<p>Renaming of function ANORM from:<br>

     Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
     Package of Special Function Routines and Test Drivers"
     acm Transactions on Mathematical Software. 19, 22-32.<br><br>

     with slight modifications to return ccum and to deal with
     machine constants.</p>
<p><strong>Original Comments:</strong>
<p>This function evaluates the normal distribution function:</p>
<p><pre>
 
                              / x
                     1       |       -t*t/2
          P(x) = ----------- |      e       dt
                 sqrt(2 pi)  |
                             /-oo
</pre></p>
<p>   The main computation evaluates near-minimax approximations
   derived from those in "Rational Chebyshev approximations for
   the error function" by W. J. Cody, Math. Comp., 1969, 631-637.
   This transportable program uses rational functions that
   theoretically approximate the normal distribution function to
   at least 18 significant decimal digits.  The accuracy achieved
   depends on the arithmetic system, the compiler, the intrinsic
   functions, and proper selection of the machine-dependent
   constants.</p>
<p><pre>
 Explanation of machine-dependent constants.
 
   MIN   = smallest machine representable number.
 
   EPS   = argument below which anorm(x) may be represented by
           0.5  and above which  x*x  will not underflow.
           A conservative value is the largest machine number X
           such that   1.0 + X = 1.0   to machine precision.
</pre></p>
<p>Error returns<br><br>
 
  The program returns  ANORM = 0     for  ARG .LE. XLOW.</p>
 
 
<p>Intrinsic functions required are: ABS, AINT, EXP</p>

\param arg Upper limit of integration.
\param result Cumulative normal distribution.                                      
\param ccum Compliment of Cumulative normal distribution.
\author                                        
    W. J. Cody,
          Mathematics and Computer Science Division,
          Argonne National Laboratory,
          Argonne, IL 60439 
\date Latest modification: March 15, 1992
*/
void cumnor(double *arg,double *result,double *ccum);

/*! \brief CUMulative POIsson distribution
 \details
     Returns the  probability  of  S   or  fewer events in  a   Poisson
     distribution with mean XLAM.
 \param s Upper limit of cumulation of the Poisson.
 \param xlam Mean of the Poisson distribution.
 \param cum Cumulative poisson distribution.
 \param ccum Compliment of Cumulative poisson distribution.
 \author
     Uses formula  26.4.21   of   Abramowitz and  Stegun,  Handbook  of
     Mathematical   Functions  to reduce   the   cumulative Poisson  to
     the cumulative chi-square distribution.
*/
void cumpoi(double *s,double *xlam,double *cum,double *ccum);

/*! \brief CUMulative T-distribution
\details Computes the integral from -infinity to T of the t-density.
\param t Upper limit of integration of the t-density.
\param df Degrees of freedom of the t-distribution. 
\param cum Cumulative t-distribution.
\param ccum Compliment of Cumulative t-distribution. 
\author
     Formula 26.5.27   of     Abramowitz  and   Stegun,    Handbook  of
     Mathematical Functions  is   used   to  reduce the  t-distribution
     to an incomplete beta. */
void cumt(double *t,double *df,double *cum,double *ccum);

/*!
\brief Double Precision Sterling Remainder for Complete
                    Beta Function
\details 
     <p>Log(Beta(A,B)) = Lgamma(A) + Lgamma(B) - Lgamma(A+B)
     where Lgamma is the log of the (complete) gamma function</p>
 
     <p>Let ZZ be approximation obtained if each log gamma is approximated
     by Sterling's formula, i.e.,
     Sterling(Z) = LOG( SQRT( 2*PI ) ) + ( Z-0.5 ) * LOG( Z ) - Z</p>
\return Log(Beta(A,B)) - ZZ
\param a One argument of the Beta
\param b The other argument of the Beta */
double dbetrm(double *a,double *b);

/*! \brief Double precision EVALuate a PoLynomial at X
 \details 
     returns
          A(1) + A(2)*X + ... + A(N)*X**(N-1)
 \param a  Array of coefficients of the polynomial.
                                        A is DOUBLE PRECISION(N)
 \param n Length of A, also degree of polynomial - 1.
 \param x Point at which the polynomial is to be evaluated. */
double devlpl(double a[],int *n,double *x);

/*! \brief Evaluation of the function EXP(X) - 1
\param x Argument at which exp(x)-1 desired
\details 
<p>Renaming of function rexp from code of below authors:-</p>
\author 
     DiDinato, A. R. and Morris,  A.   H.  Algorithm 708: Significant
     Digit Computation of the Incomplete  Beta  Function Ratios.  ACM
     Trans. Math.  Softw. 18 (1993), 360-373.
*/
double dexpm1(double*x);

/*! \brief Double precision NoRmal distribution INVerse
\details 
<p>Returns X  such that CUMNOR(X)  =   P,  i.e., the  integral from -
     infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU is P</p>
<p>N/B If P or Q .lt. machine EPS returns +/- DINVNR(EPS)</p>
\author 
  The  rational   function   on  page 95    of Kennedy  and  Gentle,
     Statistical Computing, Marcel Dekker, NY , 1980 is used as a start
     value for the Newton method of finding roots.
\param p The probability whose normal deviate is sought.
\param q 1-P                    
*/
double dinvnr(double *p,double *q);

/*! \brief DEFINE DINVR */
static void E0000(int,int*,double*,double*,unsigned long*,
                  unsigned long*,double*,double*,double*,
                  double*,double*,double*,double*);
/*! 
\brief bounds the zero of the function and invokes zror Reverse Communication
\details
     Bounds the    function  and  invokes  ZROR   to perform the   zero
     finding.  STINVR  must  have   been  called  before this   routine
     in order to set its parameters.
 \param status   At the beginning of a zero finding problem, STATUS
                 should be set to 0 and INVR invoked.  (The value
                 of parameters other than X will be ignored on this cal
 
                 <p>When INVR needs the function evaluated, it will set
                 STATUS to 1 and return.  The value of the function
                 should be set in FX and INVR again called without
                 changing any of its other parameters.</p>
 
                 <p>When INVR has finished without error, it will return
                 with STATUS 0.  In that case X is approximately a root
                 of F(X).</p>
 
                 <p>If INVR cannot bound the function, it returns status
                 -1 and sets QLEFT and QHI.</p>
\param x The value of X at which F(X) is to be evaluated.                         
\param fx The value of F(X) calculated when INVR returns with STATUS = 1.
\param qleft Defined only if QMFINV returns .FALSE.  In that
          case it is .TRUE. If the stepping search terminated
          unsucessfully at SMALL.  If it is .FALSE. the search
          terminated unsucessfully at BIG.
                    QLEFT is LOGICAL
\param qhi Defined only if QMFINV returns .FALSE.  In that
          case it is .TRUE. if F(X) .GT. Y at the termination
          of the search and .FALSE. if F(X) .LT. Y at the
          termination of the search.
                    QHI is LOGICAL
 */
void dinvr(int *status,double *x,double *fx, unsigned long *qleft,unsigned long *qhi);

/*! \brief SeT INverse finder - Reverse Communication Function
\details
<p>Given a monotone function F finds X
     such that F(X) = Y.  Uses Reverse communication -- see invr.
     This routine sets quantities needed by INVR.
          More Precise Description of INVR -
     F must be a monotone function, the results of QMFINV are
     otherwise undefined.  QINCR must be .TRUE. if F is non-
     decreasing and .FALSE. if F is non-increasing.
     QMFINV will return .TRUE. if and only if F(SMALL) and
     F(BIG) bracket Y, i. e.,
          QINCR is .TRUE. and F(SMALL).LE.Y.LE.F(BIG) or
          QINCR is .FALSE. and F(BIG).LE.Y.LE.F(SMALL)<br>
     if QMFINV returns .TRUE., then the X returned satisfies
     the following condition.  let
               TOL(X) = MAX(ABSTOL,RELTOL*ABS(X)) <br>
     then if QINCR is .TRUE.,<br>
          F(X-TOL(X)) .LE. Y .LE. F(X+TOL(X)) <br>
     and if QINCR is .FALSE.<br>
          F(X-TOL(X)) .GE. Y .GE. F(X+TOL(X))<br></P>
                              Arguments
    <p>
     Compares F(X) with Y for the input value of X then uses QINCR
     to determine whether to step left or right to bound the
     desired x.  the initial step size is
          MAX(ABSSTP,RELSTP*ABS(S)) for the input value of X.</P>
     <p>Iteratively steps right or left until it bounds X.
     At each step which doesn't bound X, the step size is doubled.
     The routine is careful never to step beyond SMALL or BIG.  If
     it hasn't bounded X at SMALL or BIG, QMFINV returns .FALSE.
     after setting QLEFT and QHI.</P>
     <p>If X is successfully bounded then Algorithm R of the paper
     'Two Efficient Algorithms with Guaranteed Convergence for
     Finding a Zero of a Function' by J. C. P. Bus and
     T. J. Dekker in ACM Transactions on Mathematical
     Software, Volume 1, No. 4 page 330 (DEC. '75) is employed
     to find the zero of the function F(X)-Y. This is routine
     QRZERO. </P>
\param zsmall The left endpoint of the interval to be
          searched for a solution.
\param zbig The right endpoint of the interval to be
          searched for a solution.
                    
\param zabsst The initial step size in the search
          is MAX(ABSSTP,RELSTP*ABS(X)). See algorithm.
\param zrelst The initial step size in the search
          is MAX(ABSSTP,RELSTP*ABS(X)). See algorithm.
\param zstpmu When a step doesn't bound the zero, the step
                size is multiplied by STPMUL and another step
                taken.  A popular value is 2.0
\param zabsto numbers that determine the accuracy
          of the solution. See function for a precise definition.
\param zrelto numbers that determine the accuracy
          of the solution. See function for a precise definition.
*/
void dstinv(double *zsmall,double *zbig,double *zabsst,
	    double *zrelst,double *zstpmu,double *zabsto,
	    double *zrelto);

/*! \brief Double precision Logarith of the Asymptotic Normal
 \details
 <p>
       Computes the logarithm of the cumulative normal distribution
      from abs( x ) to infinity for abs( x ) >= 5.</P>
 <p>23 term expansion of formula 26.2.12 of Abramowitz and Stegun.
      The relative error at X = 5 is about 0.5E-5.</P>
 <p>N/B ABS(X) must be >= 5 else there is an error stop.</P>
                              
\param x Value at which cumulative normal to be evaluated */
double dlanor(double *x);
 
/*! \brief Double precision LN(1-X)
\details 
     <p>Returns ln(1-x) for small x (good accuracy if x .le. 0.1).</p>
     <p>Note that the obvious code of
               LOG(1.0-X)
     won't work for small X because 1.0-X loses accuracy</P>
	 <p>        
     If X > 0.1, the obvious code above is used ELSE
     The Taylor series for 1-x is expanded to 20 terms.</p>

\param x Value for which ln(1-x) is desired.
*/
double dln1mx(double*x);

/*! \brief Double precision LN(1+X)
\details     
     Returns ln(1+x)<br>
     Note that the obvious code of
               LOG(1.0+X)
     won't work for small X because 1.0+X loses accuracy
 \param x Value for which ln(1-x) is desired.
 \author
     Renames ALNREL from:<br>
     DiDinato, A. R. and Morris,  A.   H.  Algorithm 708: Significant
     Digit Computation of the Incomplete  Beta  Function Ratios.  ACM
     Trans. Math.  Softw. 18 (1993), 360-373.
*/
double dln1px(double*x);

/*! \brief Double precision LN of the complete BETa
 \details  Returns the natural log of the complete beta function,
     i.e.,<br>
                  ln( Gamma(a)*Gamma(b) / Gamma(a+b)
\param a0 symmetric arguments to the complete beta
\param b0 symmetric arguments to the complete beta 
\author  
     Renames BETALN from:<br>
     DiDinato, A. R. and Morris,  A.   H.  Algorithm 708: Significant
     Digit Computation of the Incomplete  Beta  Function Ratios.  ACM
     Trans. Math.  Softw. 18 (1993), 360-373. 
*/
double dlnbet(double *a0,double *b0);

/*! \brief Double precision LN of the GAMma function 
\details         
     Renames GAMLN from:<br>
     DiDinato, A. R. and Morris,  A.   H.  Algorithm 708: Significant
     Digit Computation of the Incomplete  Beta  Function Ratios.  ACM
     Trans. Math.  Softw. 18 (1993), 360-373.
\return The natural logarithm of GAMMA(X).
\param a value at which scaled log gamma is to be returned
*/
double dlngam(double*a);

/*! \brief Double precision Sterling Remainder Function
\details                             
If Z >= 6 uses 9 terms of series in Bernoulli numbers
     (Values calculated using Maple) <br>
     Otherwise computes difference explicitly
\return
     Log(Gamma(Z))  -  Sterling(Z)  where   Sterling(Z)  is
     Sterling's Approximation to Log(Gamma(Z))<br>
     Sterling(Z) = LOG( SQRT( 2*PI ) ) + ( Z-0.5 ) * LOG( Z ) - Z
\param z Value at which Sterling remainder calculated.
           Must be positive.
*/
double dstrem(double*z);

/*! \brief Double precision Initalize Approximation to
           INVerse of the cumulative T distribution 
\details Returns  the  inverse   of  the T   distribution   function, i.e.,
     the integral from 0 to INVT of the T density is P. This is an
     initial approximation
\param p p-value whose inverse from the T distribution is
          desired.
\param q 1-p
\param df Degrees of freedom of the T distribution.
*/
double dt1(double *p,double *q,double *df);

/*! \brief DEFINE DZROR */
static void E0001(int IENTRY,int *status,double *x,double *fx,
		  double *xlo,double *xhi,unsigned long *qleft,
		  unsigned long *qhi,double *zabstl,double *zreltl,
		  double *zxhi,double *zxlo);

/*! \brief  Double precision ZeRo of a function -- Reverse Communication 
\details Performs the zero finding.  STZROR must have been called before
     this routine in order to set its parameters.
\param status  At the beginning of a zero finding problem, STATUS
                 should be set to 0 and ZROR invoked.  (The value
                 of other parameters will be ignored on this call.)<br><br>
 
                 When ZROR needs the function evaluated, it will set
                 STATUS to 1 and return.  The value of the function
                 should be set in FX and ZROR again called without
                 changing any of its other parameters.<br><br>
 
                 When ZROR has finished without error, it will return
                 with STATUS 0.  In that case (XLO,XHI) bound the answer<br><br>
 
                 If ZROR finds an error (which implies that F(XLO)-Y an
                 F(XHI)-Y have the same sign, it returns STATUS -1.  In
                 this case, XLO and XHI are undefined.
 \param x The value of X at which F(X) is to be evaluated.
 \param fx The value of F(X) calculated when ZROR returns with
            STATUS = 1.
\param xlo When ZROR returns with STATUS = 0, XLO bounds the
             inverval in X containing the solution below.
\param xhi When ZROR returns with STATUS = 0, XHI bounds the
             inverval in X containing the solution above.
\param qleft .TRUE. if the stepping search terminated unsucessfully at XLO.  If it is .FALSE. the search terminated
                unsucessfully at XHI.
\param qhi .TRUE. if F(X) .GT. Y at the termination of the
              search and .FALSE. if F(X) .LT. Y at the
              termination of the search.
*/
void dzror(int *status,double *x,double *fx,double *xlo,
	   double *xhi,unsigned long *qleft,unsigned long *qhi);

/*! \brief Double precision SeT ZeRo finder - Reverse communication version Function
\details Sets quantities needed by ZROR.<br>  The function of ZROR
     and the quantities set is given here.<br>
     Concise Description - Given a function F<br>
     find XLO such that F(XLO) = 0.<br>
          More Precise Description -
     Input condition. F is a double precision function of a single
     double precision argument and XLO and XHI are such that
          F(XLO)*F(XHI)  .LE.  0.0<br>
     If the input condition is met, QRZERO returns .TRUE.
     and output values of XLO and XHI satisfy the following<br>
          F(XLO)*F(XHI)  .LE. 0.<br>
          ABS(F(XLO)  .LE. ABS(F(XHI)<br>
          ABS(XLO-XHI)  .LE. TOL(X)<br>
     where<br>
          TOL(X) = MAX(ABSTOL,RELTOL*ABS(X))<br><br>
     If this algorithm does not find XLO and XHI satisfying
     these conditions then QRZERO returns .FALSE.  This
     implies that the input condition was not met.<br>
 \param zxlo  The left endpoint of the interval to be
           searched for a solution.
                    XLO is DOUBLE PRECISION
\param zxhi The right endpoint of the interval to be
           for a solution.
                    XHI is DOUBLE PRECISION
\param zabstl etermine the accuracy
                      of the solution. See function for a
                      precise definition.
\param zreltl etermine the accuracy
                      of the solution. See function for a
                      precise definition.
\author
Algorithm R of the paper 'Two Efficient Algorithms with
     Guaranteed Convergence for Finding a Zero of a Function'
     by J. C. P. Bus and T. J. Dekker in ACM Transactions on
     Mathematical Software, Volume 1, no. 4 page 330
     (Dec. '75) is employed to find the zero of F(X)-Y.
*/
void dstzr(double *zxlo,double *zxhi,double *zabstl,double *zreltl);

/*! \brief EVALUATION OF THE double ERROR FUNCTION */
double erf1(double*);

/*! \brief EVALUATION OF THE COMPLEMENTARY ERROR FUNCTION
\details  ERFC1(IND,X) = ERFC(X)            IF IND = 0<br>
          ERFC1(IND,X) = EXP(X*X)*ERFC(X)   OTHERWISE
*/
double erfc1(int *ind,double *x);

/*! \brief EVALUATION OF EXP(MU + X) */
double esum(int *mu,double *x);

/*! \details <p>IF L = 0 THEN  EXPARG(L) = THE LARGEST POSITIVE W FOR WHICH
     EXP(W) CAN BE COMPUTED.</P>
 
     <p>IF L IS NONZERO THEN  EXPARG(L) = THE LARGEST NEGATIVE W FOR
     WHICH THE COMPUTED VALUE OF EXP(W) IS NONZERO.</p>
 
     <p>N/B ... ONLY AN APPROXIMATE VALUE FOR EXPARG(L) IS NEEDED.</P>*/
double exparg(int*l);

/*! \brief EVALUATION OF I (A,B) X FOR B .LT. MIN(EPS,EPS*A) AND X .LE. 0.5.
\details SET  FPSER = X**A
*/ 
double fpser(double *a,double *b,double *x,double *eps);

/*! \brief COMPUTATION OF 1/GAMMA(A+1) - 1  FOR -0.5 .LE. A .LE. 1.5 */
double gam1(double*a);

/*! \brief INVERSE INCOMPLETE GAMMA RATIO FUNCTION
\details 
     <p>GIVEN POSITIVE A, AND NONEGATIVE P AND Q WHERE P + Q = 1.
     THEN X IS COMPUTED WHERE P(A,X) = P AND Q(A,X) = Q. SCHRODER
     ITERATION IS EMPLOYED. THE ROUTINE ATTEMPTS TO COMPUTE X
     TO 10 SIGNIFICANT DIGITS IF THIS IS POSSIBLE FOR THE
     PARTICULAR COMPUTER ARITHMETIC BEING USED.</p>
 
                      
 
     <p>X IS A VARIABLE. IF P = 0 THEN X IS ASSIGNED THE VALUE 0,
     AND IF Q = 0 THEN X IS SET TO THE LARGEST FLOATING POINT
     NUMBER AVAILABLE. OTHERWISE, GAMINV ATTEMPTS TO OBTAIN
     A SOLUTION FOR P(A,X) = P AND Q(A,X) = Q. IF THE ROUTINE
     IS SUCCESSFUL THEN THE SOLUTION IS STORED IN X.</P>
 
     <p>X0 IS AN OPTIONAL INITIAL APPROXIMATION FOR X. IF THE USER
     DOES NOT WISH TO SUPPLY AN INITIAL APPROXIMATION, THEN SET
     X0 .LE. 0.</p>
 
     <p>IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
     WHEN THE ROUTINE TERMINATES, IERR HAS ONE OF THE FOLLOWING
     VALUES ...</p>
 <p>
       IERR =  0    THE SOLUTION WAS OBTAINED. ITERATION WAS
                    NOT USED.<br>
       IERR.GT.0    THE SOLUTION WAS OBTAINED. IERR ITERATIONS
                    WERE PERFORMED<br>.
       IERR = -2    (INPUT ERROR) A .LE. 0<br>
       IERR = -3    NO SOLUTION WAS OBTAINED. THE RATIO Q/A
                    IS TOO LARGE.<br>
       IERR = -4    (INPUT ERROR) P + Q .NE. 1<br>
       IERR = -6    20 ITERATIONS WERE PERFORMED. THE MOST
                    RECENT VALUE OBTAINED FOR X IS GIVEN.
                    THIS CANNOT OCCUR IF X0 .LE. 0.<br>
       IERR = -7    ITERATION FAILED. NO VALUE IS GIVEN FOR X.
                    THIS MAY OCCUR WHEN X IS APPROXIMATELY 0.<br>
       IERR = -8    A VALUE FOR X HAS BEEN OBTAINED, BUT THE
                    ROUTINE IS NOT CERTAIN OF ITS ACCURACY.
                    ITERATION CANNOT BE PERFORMED IN THIS
                    CASE. IF X0 .LE. 0, THIS CAN OCCUR ONLY
                    WHEN P OR Q IS APPROXIMATELY 0. IF X0 IS
                    POSITIVE THEN THIS CAN OCCUR WHEN A IS
                    EXCEEDINGLY CLOSE TO X AND A IS EXTREMELY
                    LARGE (SAY A .GE. 1.E20).
					</p>
 \author
     ALFRED H. MORRIS, JR.
        NAVAL SURFACE WEAPONS CENTER
        DAHLGREN, VIRGINIA */
void gaminv(double *a,double *x,double *x0,double *p,double *q, int *ierr);

/*! \brief EVALUATION OF LN(GAMMA(A)) FOR POSITIVE A
\details D = 0.5*(LN(2*PI) - 1)
\author ALFRED H. MORRIS
          NAVAL SURFACE WARFARE CENTER
          DAHLGREN, VIRGINIA */
double gamln(double*);

/*! \brief EVALUATION OF LN(GAMMA(1 + A)) FOR -0.2 .LE. A .LE. 1.25 */
double gamln1(double *a);

/*! \brief EVALUATION OF THE GAMMA FUNCTION FOR double ARGUMENTS
\details GAMMA(A) IS ASSIGNED THE VALUE 0 WHEN THE GAMMA FUNCTION CANNOT
     BE COMPUTED.
\author ALFRED H. MORRIS, JR.
          NAVAL SURFACE WEAPONS CENTER
          DAHLGREN, VIRGINIA
*/
double Xgamm(double*a);

/*! \brief EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS
                      P(A,X) AND Q(A,X)
\details
     IT IS ASSUMED THAT A .LE. 1.  EPS IS THE TOLERANCE TO BE USED.
     THE INPUT ARGUMENT R HAS THE VALUE E**(-X)*X**A/GAMMA(A). */
void grat1(double *a,double *x,double *r,double *p,double *q, double *eps);

/*! \brief EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS
                      P(A,X) AND Q(A,X)
\details <p>IT IS ASSUMED THAT A AND X ARE NONNEGATIVE, WHERE A AND X
     ARE NOT BOTH 0.</p>
 
     <p>ANS AND QANS ARE VARIABLES. GRATIO ASSIGNS ANS THE VALUE
     P(A,X) AND QANS THE VALUE Q(A,X). IND MAY BE ANY INTEGER.
     IF IND = 0 THEN THE USER IS REQUESTING AS MUCH ACCURACY AS
     POSSIBLE (UP TO 14 SIGNIFICANT DIGITS). OTHERWISE, IF
     IND = 1 THEN ACCURACY IS REQUESTED TO WITHIN 1 UNIT OF THE
     6-TH SIGNIFICANT DIGIT, AND IF IND .NE. 0,1 THEN ACCURACY
     IS REQUESTED TO WITHIN 1 UNIT OF THE 3RD SIGNIFICANT DIGIT.</p>
 
     <p>ERROR RETURN ...
        ANS IS ASSIGNED THE VALUE 2 WHEN A OR X IS NEGATIVE,
     WHEN A*X = 0, OR WHEN P(A,X) AND Q(A,X) ARE INDETERMINANT.
     P(A,X) AND Q(A,X) ARE COMPUTATIONALLY INDETERMINANT WHEN
     X IS EXCEEDINGLY CLOSE TO A AND A IS EXTREMELY LARGE.</p>
 \author ALFRED H. MORRIS, JR.
        NAVAL SURFACE WEAPONS CENTER
        DAHLGREN, VIRGINIA
 */
void gratio(double *a,double *x,double *ans,double *qans,int *ind);

/*! \brief EVALUATION OF THE FUNCTION LN(GAMMA(A + B))
          FOR 1 .LE. A .LE. 2  AND  1 .LE. B .LE. 2 */
double gsumln(double*a,double*b);

/*! \brief EVALUATION OF THE DIGAMMA FUNCTION
 \details
     <p>PSI(XX) IS ASSIGNED THE VALUE 0 WHEN THE DIGAMMA FUNCTION CANNOT
     BE COMPUTED.</P>
 
     <p>THE MAIN COMPUTATION INVOLVES EVALUATION OF RATIONAL CHEBYSHEV
     APPROXIMATIONS PUBLISHED IN MATH. COMP. 27, 123-127(1973) BY
     CODY, STRECOK AND THACHER.</p>
 
\authors
     PSI WAS WRITTEN AT ARGONNE NATIONAL LABORATORY FOR THE FUNPACK
     PACKAGE OF SPECIAL FUNCTION SUBROUTINES. PSI WAS MODIFIED BY
     A.H. MORRIS (NSWC).*/
double psi(double *xx);

/*! \brief EVALUATION OF EXP(-X)*X**A/GAMMA(A)
\details RT2PIN = 1/SQRT(2*PI) */
double rcomp(double*a,double*x);

/*! \brief EVALUATION OF THE FUNCTION EXP(X) - 1 */
double rexp(double*x);

/*! \brief COMPUTATION OF  X - 1 - LN(X) */
double rlog(double*x);

/*! \brief EVALUATION OF THE FUNCTION X - LN(1 + X) */
double rlog1(double *x);

/*!
\brief PROVIDES SINGLE PRECISION MACHINE CONSTANTS
\details 
     <p>SPMPAR PROVIDES THE SINGLE PRECISION MACHINE CONSTANTS FOR
     THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT
     I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE
     SINGLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND
     ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN</P>
 
        <p>SPMPAR(1) = B**(1 - M), THE MACHINE PRECISION,<br>
 
        SPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE,<br>
 
        SPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE.<br></p>
 \author ALFRED H. MORRIS, JR. (NAVAL SURFACE WARFARE CENTER, DAHLGREN, VIRGINIA)<br>
		 MODIFIED BY BARRY W. BROWN */
double spmpar(int *i);

/*!
\brief STarting VALue for Newton-Raphon calculation of Normal distribution Inverse
\details <p>Returns X  such that CUMNOR(X)  =   P,  i.e., the  integral from -
     infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU is P</p>
	 <p>The  rational   function   on  page 95    of Kennedy  and  Gentle,
     Statistical Computing, Marcel Dekker, NY , 1980.</p>
 \param p The probability whose normal deviate is sought. P is DOUBLE PRECISION */
 double stvaln(double *p);

/*! \brief Truncates a double precision number to an integer and returns the value in a double. */
double fifdint(double a);

/*! \brief returns the maximum of two numbers a and b */
double fifdmax1(double a,double b);

/*! \brief returns the minimum of two numbers a and b */
double fifdmin1(double a,double b);

/*! \brief transfers the sign of the variable "sign" to the variable "mag" */
double fifdsign(double mag,double sign);

/*! \brief Truncates a double precision number to a long integer */
long fifidint(double a);

/*! \brief returns the modulo of a and b */
long fifmod(long a,long b);

/*! \brief Prints msg to standard error and then exits */
void ftnstop(char* msg);

#ifdef __cplusplus
}
#endif

#endif // _CDFLIB_H_
